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ABSTRACT 

Due to the absence of direct measurement, the magnetic field in the solar corona is usually ex- 
trapolated from the photosphere in numerical way. At the moment, the nonlinear force- free field 
(NLFFF) model dominates the physical models for field extrapolation in the low corona. Recently 
we have developed a new NLFFF model with MHD relaxation to reconstruct the coronal magnetic 
field. This method is based on CESE-MHD model with the conservation-element/solution-element 
(CESE) spacetimc scheme. In this paper, we report the application of the CESE-MHD-NLFFF 
code to 5'£'0/IIMI data with magnetograms sampled for two active regions (ARs), NOAA AR 11158 
and 11283, both of which were very non-potential, producing X-class flares and eruptions. The raw 
magnetograms are preprocessed to remove the force and then inputted into the extrapolation code. 
Qualitative comparison of the results with the SDO/AIA images shows that our code can reconstruct 
magnetic field lines resembling the EUV-observed coronal loops. Most important structures of the 
active regions are reproduced excellently, like the highly-sheared field lines that suspend filaments in 
AR 11158 and twisted flux rope which corresponds to a sigmoid in AR 11283. Quantitative assess 
of the results shows that the force-free constraint is fulfilled very well in the strong-field regions but 
apparently not that well in the weak-field regions because of data noise and numerical errors in the 
small currents. 
Subject headings: Magnetic fields; Magnetohydrodynamics (MHD); Methods: numerical; Sun: corona 



1. INTRODUCTION 

Magnetic field extrapolation is an important tool to 
study the three-dimensional (3D) solar coronal magnetic 
field, which is still difficult to measure directly (Saku- 
rai 1989; Aly 1989; Amari et al. 1997; McClymont et al. 
1997; Wiegelmann 2008; DeRosa et al. 2009). The mod- 
els being used most popularly for field extrapolation 
are the potential field model, the linear force-free field 
model, and the nonlinear force-free field (NLFFF) model. 
These models are all based on the same assumption that 
the Lorentz force is self-balancing in the corona, but 
adopt different simplifications of the current distribu- 
tion. Among these models. The NLFFF is the most used 
one for characterizing magnetic field in the low corona, 
where there is significant and localized electric current, 
especially in the active regions (ARs). 

But, to directly solve the general NLFFF equation 



(V X B) X B = 0, V • B = 



(1) 



is really difficult. As is known, the system is nonlinear 
intrinsically and even the existence and uniqueness of a 
solution for a given boundary condition are not proved 
theoretically; solutions have rarely been found in closed 
analytic form (e.g.. Low & Lou 1990) and in most cases 
people can only resort to numerical method using com- 
puter (many numerical codes have been developed in the 
past decades, e.g., Wu ct al. 1990; Roumcliotis 1996; 
Amari et al. 1999; Wheatland et al. 2000; Yan & Sakurai 
2000; Wiegelmann 2004; Valori et al. 2007; Inoue et al. 
2011, one may refer to a recent living review by Wiegel- 
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mann & Sakurai (2012)). Moreover the observation can 
only provide a bottom boundary of data, and even worse, 
on the photosphere the field is forced significantly by the 
dense plasma and thus conflicts with the fundamental 
force-free assumption. Besides the noise in the obser- 
vation, measurement error and instrumental uncertainty 
(e.g., the well known 180° ambiguity of the transverse 
fields) are all rather unfavorable for practical computa- 
tion. Thus the observed magnetogram usually needs to 
be preprocessed to remove the force and noise for provid- 
ing a better input (Wiegelmann et al. 2006). Anyhow, at 
present one can hardly seek an exact force-free solution 
with the observation information fully satisfied. The best 
we can do is to find a good balance between the force-free 
constraint and deviation from the real observation, i.e., 
to seek an approximately force-free solution that matches 
the photospheric field measurements as well as possible. 
Recently we have developed a new extrapolation code 
called CESE-MHD-NLFFF- (Jiang et al. 2011; Jiang & 
Feng 2012a), which is based on magnetohydrodynamics 
(MHD) relaxation method and an advanced numerical 
scheme, the spacetime conservation-element/solution- 
element (CESE) method, for faster convergence and bet- 
ter accuracy over the available codes. The good perfor- 
mance and high accuracy of the code have been demon- 
strated through critical comparisons with previous joint 
studies by Schrijver et al. (2006) and Metcalf et al. 

^ We initially planned to develop a full MHD model for comput- 
ing both the static non-potential field and dynamic evolution of 
ARs (i.e., the CESE-MHD model, Jiang et al. 2011). But wc find 
that the MHD solver is rather slow to construct the static field, 
although success has been reported in Jiang et al. (2012a). We 
thus move to develop a NLFFF version of this model for a faster 
convergence speed, while the full MHD solver is more suitable for 
simulating transient events like the eruptions. 
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(2008), in which various NLFFF codes are assessed based 
on several NLFFF benchmark tests. We have also suc- 
cessfully extended the code to application in spherical 
geometry and seamless full-sphere extrapolation for the 
global corona (Jiang ct al. 2012b). 

In this paper we report the application of the CESE- 
MHD-NLFFF code to real solar data, i.e., the presently 
released SDO/BMl magnetograms. To deal with the real 
observation data, we also have developed a new prepro- 
cessing method to remove the force in the raw magne- 
togram (Jiang & Feng 2013). Extra advancements are 
made to the original code to further enhance the ability 
of handling the high-resolution but noisy data. Magne- 
tograms of two ARs, AR 11158 and AR 11283 are sam- 
pled for our tests of extrapolation. The results are care- 
fully assessed both qualitatively and quantitatively. We 
show that our code can recover magnetic field lines re- 
sembling the plasma loops seen in the SDO/AIA images, 
and reproduce most important structures of the ARs re- 
markably well like the highly-sheared field lines that sus- 
pend filaments in AR 11158 and twisted flux rope that 
corresponds to a sigmoid in AR 11283. 

The remainder of the paper is organized as follows. We 
first describe briefly the CESE-MHD-NLFFF code along 
with its improvements in Section 2. The magnetogram 
from SDO/HMI and the preprocessed result of the raw 
data are given in Section 3. We then present the ex- 
trapolation results for these data including both the raw 
and preprocessed magnetogram in Section 4. Finally, we 
draw conclusions and give some outlooks for future work 
in Section 5. 

2. THE CESE-MHD-NLFFF CODE 

The basic idea of using the MHD relaxation approach 
to solve the force-free field is to use some kind of ficti- 
tious dissipation to drive the MHD system to an equi- 
librium in which all the forces can be neglected compar- 
ing with the Lorentz force while the boundary magne- 
togram is satisfied. In this way the Lorentz force should 
be self-balancing and the field can be regarded as the tar- 
get force-free solution. We solve the magneto-frictional 
model equations in the magnetic splitting form as 



dpv 

1h 



= (V X Bi) X B - (V • Bi)B - i/pv, 



aBi 



= V X (v X B) + V(/^V • Bi) - vV • Bi, 

p=|B|2 + po, B = Bo + Bi. (2) 

B is the target force-free field to be solved, Bq is the 
potential field matching the normal component of the 
magnetogram, and Bi is the deviation between B and 
Bq. V is the frictional coefficient and /i is the numerical 
diffusive speed of the magnetic monopole. The value of 
them are respectively given by i^ = l/(5Ai) and /i = 
0.4(Ax)^/At in the code, according to the time step Ai 
and local grid size Az. Many advantages can be gained 
by solving such form of above equations (Jiang & Feng 
2012a; Jiang et al. 2012b). 

The above equation system (2) is solved by our CESE- 
MHD scheme (Jiang ct al. 2010). In principle we can use 
any available MHD code to solve this set of equations, 
since it is a subset of the full MHD system. Taking con- 
sideration of the computational efficiency and accuracy. 



we prefer to utilize modern advanced MHD codes. How- 
ever, most of the modern MHD codes are based on theory 
of characteristic decomposition of a hyperbolic system, 
thus are not suitable for equation (2), which is not a 
hyperbolic system. The CESE scheme is a new method 
free of characteristic decomposition and is very suitable 
for the equation form here. Furthermore, the CESE- 
MHD code has been extensively used in solar physics, 
e.g., the data-driven evolution modeling of AR (Jiang 
et al. 2012a), the global corona (Feng et al. 2012) and 
the interplanetary solar wind (Feng et al. 2012; Yang 
et al. 2012). 

To adapt for the application to real solar data, we 
have made extra improvements to the previous version 
of CESE-MHD-NLFFF. The first improvement is made 
to enhance the ability of handling noisy data in the mag- 
netograms. In the noisy weak-field regions of magne- 
tograms (where the signal-to-noise ratio is small, say 
|B| < 100 G), the term (V x B) x B/|B|2 could be very 
large due to numerical gradients of the random noise, 
thus the velocity v is prone to be accelerated to extremely 
high, which can severely restrict the time step and slow 
the relaxation process of the entire system, even mak- 
ing the computation unmanageable. To deal with this 
difficulty, the pseudo-plasma density p is designed with 
pQ — i?^^;jjexp(— z/iJjji), where z is the height from the 
bottom surface and -Bmin = 100 G, 7?„i = 5 pixel. In this 
way, the velocity (near the bottom magnetogram) in the 
weak-field regions can be reduced significantly, while it 
is barely affected in the strong-field regions. 




Fig. 1. — The grid structure: the entire volume is divided into 
blocks and each block has 8x8x8 cells. Slices through the vol- 
ume in three axis directions are plotted to show the structure of 
the blocks and the bottom contour map represents B^ on the pho- 
tosphere. 



Secondly, to deal with the high-resolution observation 
data, the extrapolation is performed on a non-uniform 
grid within a block-structured, distributed-memory par- 
allel computational framework (e.g., Jiang et al. 2012a). 
Specifically, the whole computational volume is divided 
into blocks with different spatial resolution, and the 
blocks are evenly distributed among the processors. 
Within this framework, we have lots of freedom to con- 
figure the mesh and save computational resources com- 
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paring with a uniform grid. As concentrated strongly in 
the photosphere but expanding rapidly into the corona 
due to abrupt drop of the gas pressure, the coronal field 
becomes smoother and weaker successively with height. 
Naturally we use a grid with decreasing resolution with 
height: at the bottom the grid spacing matches the reso- 
lution of the magnetogram and up to the top of the model 
box, the grid spacing is increased by, say, four times. 
At present we use the same resolution in the horizontal 
plane, and application of adaptive resolution based on 
the pattern of magnetic flux distribution is under devel- 
opment. With this framework, we also add some coarse 
buffer blocks around the central volume to reduce influ- 
ence by the numerical boundaries without adding much 
computational burden. An example of the grid is shown 
in Figure 1 . 

Routinely the quantities in the computational volume 
is initialized by setting Bi = and v = 0. The poten- 
tial field Bq is obtained by a Green's function method 
(e.g., Metcalf et al. 2008). The bottom boundary is in- 
crementally fed with the observed vector magnetogram 
in tens of Alfven time ta , while all the other numerical 
boundaries are fixed with Bi = and v = 0. 

3. DATA 
3.1. The HMI Data 

The Helioseismic and Magnetic Imager (HMI) on 
board the Solar Dynamics Observatory (SDO) provides 
photospheric vector magnetograms with a high resolu- 
tion both in space and time. It observes the full Sun 
with a 4Kx4K CCD whose spatial sampling is 0.5 arcsec 
per pixel. Raw filtergrams are obtained at 6 different 
wavelengths and 6 polarization states in the Fe i 6173 A 
absorption line, and are collected and converted to ob- 
servable quantities (like Dopplergrams, continuum filter- 
grams, and line-of-sight and vector magnetograms) on a 
rapid time cadence. For the vector magnetic data, each 
set of filtergrams takes 135 s to be completed. To obtain 
vector magnetograms, Stokes parameters are first derived 
from filtergrams observed over a 12-min interval and then 
inverted through the Very Fast Inversion of the Stokes 
Vector (Borrero et al. 2011). The 180° azimuthal ambi- 
guity in the transverse field is resolved by an improved 
version of the "minimum energy" algorithm (Leka et al. 
2009). Regions of interest with strong magnetic field are 
automatically identified near real time (Turmon et al. 
2010). A detailed description on how the vector mag- 
netograms are produced can be found on the website // 
http : // j soc . Stanford . edu/ j socwiki/VectorPaper. 

The magnetogram data we use here is down- 
loaded from http://jsoc.stanford.edu/jsocwiki/ 
VectorPaper, where the HMI vector magnetic field data 
series hmi .B_720s_el5wl332 are released for several ARs. 
There are two special formats, i.e., direct cutouts and 
remapped images. We use the remapped format which 
is more suitable for modeling in local Cartesian coordi- 
nates, since the images are computed with a Lambert 
cylindrical equal area projection centered on the tracked 
region. For our test, we select two ARs, AR 11158 and 
AR 11283, both of which produced X-class fiares and 
were very non-potential. The full resolution of the data 
is about 0.5" per pixel and we rebin them to 1" per pixel 
for the NLFFF modehng. 



Bz(G) 



-1000 -500 




Fig. 2.— Vector magnetograms for AR 11158 and AR 11283. 
The background shows the vertical components with saturation 
values of ±1000 G; the vectors represent the transverse fields (above 
200 G). The length unit is arcsec. 



AR 11158 is a well-known target studied in many re- 
cent works for different purposes (e.g., Schrijver et al. 
2011; Sun et al. 2012b; Liu et al. 2012; Jing et al. 2012), 
and was selected by Wiegelmann et al. (2012) for special 
test on optimizing their extrapolation code with HMI 
data. This AR has a multipolar and complex struc- 
ture. It produced several major flares and coronal mass 
ejections (CMEs) during its disk passage, including the 
first X-class fiare of cycle 24 on 2012 February 15. Fig- 
ure 2 (the first panel) shows a vector magnetogram of AR 
11158 taken at 20:36 UT on 2011 February 14, which wiU 
be used for our computation. This AR is well isolated 
from other ones and is almost flux-balanced (see Table 1), 
with the main polarities concentrated in the central field 
of view (FoV). As can be seen, the field shows a strong 
shearing along the polarity inversion lines (PILs). 

Another one, AR 11283, is also very eruptive, generat- 
ing several X-class flares and CMEs. We select a magne- 
togram taken at 22:00 UT on September 6, just prior to 
a major flare at 22:20 UT. The magnetogram is shown in 
the second panel of Figure 2. As an input for extrapola- 
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tion, this data is not good as AR 11158's, since the flux is 
dispersed with some strong polarities almost on the edge 
of the FoV. Also the total flux is not well balanced (see 
Table 1), which is unfavorable for our extrapolation. 

3.2. Data Preprocessing 

Generally, the raw magnetogram cannot be inputted 
directly into the NLFFF code because the intrinsic non- 
force-freeness of the photosphcric field violates the force- 
free assumption. According to the derivation by Molo- 
denskii (1969) and Aly (1989), an ideally force-free mag- 
netogram should fulfill the following conditions: 



B, dx dy = 0, F,^ B,B, dx dy = 0, 

s Js 

Fy^ ByB, dx dy = 0, F,, ^ Eb dx dy ^ 0, 
JS Js 

Tx ^ yEs dx dy = 0, Ty = xEb dx dy = 0, 
Js Js 

T,= [ {yBxB, - xByB,) dx dy = 0(3) 
Js 

where Eb ^ B^ + By — B^. These expressions are de- 
rived from the volume integrals of the total divergence, 
magnetic force and torque of an ideally force-free field 

/ V-Bdy = 0, /jxBdV" = 0, /rx(jxB)dy = 
Jv Jv Jv 

. (4) 
(by using Gauss' divergence theorem, the volume inte- 
grals 4 can be transformed to the surface integrals 3). 
In the case of the coronal field, the surface integrals of 
Equation (3) are usually restricted within the bottom 
magnetogram since the contribution from other bound- 
aries can be neglected. To assess the real data with re- 
spect to the force-free condition, three parameters are 
usually computed as (Wiegelmann et al. 2006) 



Cfiu 



J„ Bz dx dy 
Js\B.\dxdy' 



efo 



l^.l 



IF, 



Jg Pb dx dy 



\n 



IT. 



^torque 



Jg \/x2 -I- y2pg da; dy 



(5) 



where Pb = B^ + By+B^. Small values of these quanti- 
ties, e.g., eflux, eforce, etorque < 1 indicate a good input for 
the NLFFF modeling. Table 1 shows that for AR 11158, 
the force-free condition is satisfied quiet well with flux 
almost balanced and eforce, etorque less than 0.1; while for 
AR 11283, it is worse. Note that the flux non-balance 
will pose a negative effect on the extrapolation (see Sec- 
tion 4). 

Besides the non-force-freeness, the observed data con- 
tains measurement noise which is also unfavorable for 
practical implementation of extrapolation. To this end, 
preprocessing of the raw magnetogram has been pro- 
posed by Wiegelmann et al. (2006) to remove the force 
and noise for providing better input for NLFFF model- 
ing. To be consistent with our extrapolation code us- 
ing a magnetic splitting form, we recently developed a 
new code for magnetogram preprocessing (Jiang & Feng 
2013), in which the vector magnetogram is also split into 
a potential field part and a non-potential field part and 



we deal with the two parts separately. Preprocessing of 
the potential part is simply performed by taking the data 
sliced at a plane about 400 km above the photosphere'^ 
from the 3D potential field which is extrapolated from 
the observed vertical field. Then the non-potential part 
is modified and smoothed by an optimization method 
to fulfill the constraints of total magnetic force-freeness 
and torque-freeness, which is similar to the method pro- 
posed by Wiegelmann et al. (2006). We have paid par- 
ticular attention to the extents the force is needed to 
be removed and the smoothing can be performed. As for 
practical computation based on numerical discretization, 
an accurate satisfaction of force-free constraints is appar- 
ently not necessary. Also the extent of the smoothing for 
the data needs to be carefully determined, if we want to 
mimic the expansion of the magnetic field from the pho- 
tosphere to some specific height above. We use the values 
of force-freeness and smoothness calculated from the pre- 
processed potential-field part as a reference to guide the 
preprocessing of the non-potential field part, i.e., we re- 
quire that the target magnetogram has the same level 
of force-freeness and smoothness as its potential part. 
These requirements can restrict well the free parameters, 
i.e., the weighted factors in the optimization function. 

The results of preprocessed data are also given in Ta- 
ble 1 together with results of their numerical potential- 
field part. For both magnetograms, the preprocessing re- 
duces the parameters eforce and etorque by more than one 
order of magnitude, making the residual force around the 
level of numerical error (i.e., the parameters are close to 
those of the numerical potential field). The parameters 
SxjSy^Sz in the table measure the smoothness of the 
components B^ [m = x, y, z), which are defined by (see 
Jiang & Feng 2013) 



5„, = ^[(AS„0']/E[(^^™)'] 



(6) 



where the summation ^ is over all the pixels of the 
magnetogram, and A is a usual five-point 2D-Laplace 
operator, i.e., for pixel {i,j) 

ABij = ^i+i,i + -Sj-ij" + ^ij+i + ^hj-i + 4i?ij-(7) 

As shown the smoothness of the preprocessed data is 
very close to those of their potential part, and is consis- 
tent among three components, which is unlike the raw 
data with very different smoothness for different com- 
ponents. Smoothing of the data can be clearly seen by 
comparing the raw and preprocessed magnetograms as 
shown in Figure 3, and especially in the J, map which 
shows the random noise is suppressed effectively (but not 
totally) . 

Readers should be reminded that the constraints of 
Equation (3) are only necessary conditions, but not suf- 
ficient for an ideally force-free magnetogram, meaning 
that the magnetogram may still contains force even these 
conditions are satisfied. What we can say is that the 

^ We choose such height because the field becomes force-free in 
the chromosphere roughly 400 km above the photosphere accord- 
ing to (Mctcalf ct al. 1995). Our preprocessing code is designed to 
modify the photosphcric magnetogram to mimic a force-free chro- 
mospheric magnetogram at such height. 
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3. — Raw and preprossessed magnctograms for AR 11158. 
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Fig. 4. — Raw and preprossessed magnetograms for AR 11283. 

preprocessed magnetogram is more suitable for NLFFF 
modeling than the raw data, but not a completely con- 
sistent boundary condition. 



4. RESULTS 
4.1. AR 11158 
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The observed coronal loops in X-ray and EUV images 
give us a proxy of the magnetic field line geometry and 
are thus a good constraint for the magnetic field model 
besides the vector magnetograms (e.g., Aschwanden et al. 
2012; Malanushenko et al. 2012). In Figure 5 we com- 
pare the extrapolated field with coronal loops observed 
by SDO/klk in the wavelength of 171 A, in which the 
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TABLE 2 

Results of the metrics for AR 11158. Full region is extrapolation box of [0, 300](2;) x [0, 300](j/) x [0, 150] (z). Region A is 

[53,222] X [104,191] x [0,50]; Region B is [116,181] x [145,175] x [0,30]; Region C is [62,102] x [132,161] x [0,20]. The FoVs of the 

regions are shown by the boxes in Figure 6. The energy unit is 10^^ erg. 



loops are visible the best than in other channels. The 
field lines are traced from the photosphere at locations 
selected roughly according to the footpoints of the vis- 
ible bright loops, and are rendered with different colors 
for a better looking. The view angle of the field lines 
is aligned with the AIA image. We plot the field lines 
and the AIA image both alone and overlaid for a better 
inspecting of the result. As shown from an overview of 
the images, most of the extrapolated field lines closely 
resemble the observed loops. At the central region the 
field lines are strongly sheared along the major PIL and 
slightly twisted (as compared with the potential solu- 
tion), indicating the existence of strong electric currents 
along the field lines. In the left panel of Figure 6, we plot 
an image of vertical integral of the current density in the 
extrapolation volume, and the strong-current regions are 
outlined by the boxes (labeled as A, B and C). Note that 
the currents are strongly localized within the central re- 
gion B and a smaller region C. The magnetic structures 
of the strong-current regions are shown in the right panel 
of Figure 6. As can be seen, the twist of the field lines 
in region C is much stronger than that in region B. The 
results support observation studies which show that a 
filament related with a X-class flare and CIVIE exists in 
the core region B (Sun et al. 2012b), and there are small 
eruptions in region C due to flux emergence (Sun et al. 
2012a). 

Nevertheless, we note that the misalignments between 
the modeling field lines and the observation are also ob- 
vious, especially, the large loops near the north-west 
boundary of the FoV. There are several reasons for the 
misalignments: first, a local Cartesian coordinate sys- 
tem is not adequate to include the large loops which 
obviously need spherical geometry; second, the FoV of 
magnetogram may be not large enough to properly char- 



acterize the entire relevant current system, which again 
needs the curvature of the Sun's surface to be taken into 
account; third, it is not easy to locate precisely the pho- 
tosphcric footpoints of different loops that spread apart 
distinctly in the corona but are rooted very closely in the 
photosphere; fourth, the coronal field may be rather dy- 
namic, e.g., expands or oscillates due to eruptions, which 
makes the static extrapolation fail. We know that this 
visual comparison between the model result and obser- 
vation is very preliminary, and more critical comparison 
is requred, for example, with the 3D loops reconstructed 
with multi-points observation (DcRosa et al. 2009; As- 
chwanden 2011). 

Routinely, we check the quality of the numerical re- 
sult by computing several metrics. The force-freeness of 
the extrapolation data is usually measured by a current- 
weighted sine metric (CWsin) defined by 



CWsin ; 






Jx B 



(8) 



where B — |B|, J — |J| and V is the computational 
volume. We also compute a current-square-weighted sine 
metric (C^Wsin) similarly defined by 



„ L .PadV 

with more weight on the strong-current regions, 
divergence-freeness is measured by {\fi\) 
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We care about different energy contents, i.e., the total 
energy -Etot, the potential energy Spot and the free en- 
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Fig. 5. — Comparison of extrapolation field lines with AIA 171 A loops for AR 11158: the NLFFF lines (a), the potential field lines (b), 
the AIA image (c) and NLFFF lines overlaying the AIA image (d). Contour lines for ±1000 G (the black curves) of LoS photospheric field 
are over-plotted on the AIA images, and for all the panels the field lines are traced from the same set of footpoints on the bottom surface. 
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(11) 
where i?pot is the potential field strength. Results of the 
metrics are given in Table 2 for extrapolations from both 
the raw and preprocessed magnetograms. We compute 
the metrics for four different regions including the full 
extrapolation box and the subregions A, B, and C as 
outlined in Figure 6. 



For the full region, our results of the current-weighted 
sine is ~ 0.3, which means the mean misalignment an- 
gle between the magnetic field and current is about 17°. 
Such value is much larger than those from our previous 
benchmark tests using ideal or synthetic magnetograms 
(which are ^ 0.1 (6°) or smaller, see Jiang & Feng 2012a), 
but is comparable to previously reported results by other 
NLFFF codes on real magnetograms (e.g., the average 
CWsin by various NLFFF codes apphed to AR 10930 
(Schrijver et al. 2008) and AR 10953 (DeRosa et al. 2009) 
are 0.36 and 0.28, respectively). With such a large mis- 
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100 150 200 

X (arcsec) 

Fig. 6. — Strong-current regions and their magnetic structures for AR 11158. Left panel is image of the vertical integral of current density, 
i.e., Jc = J \J\dz where the current is calculated by J = V X B with a unit of G pixel"^; Regions of strong current are outlines by the 
boxes A, B and C. Right panels show the magnetic structures of the strong-current regions of B and C. 



alignment angle, the result seems to be far away from 
an exactly force-free solution which has a zero misalign- 
ment angle. However it should be noted that the metric 
CWsin may not be a good monitor for numerical so- 
lutions, which unavoidably have random numerical er- 
rors because of limited resolution. As a simple example, 
CWsin is close to 1 even for a potential field solution 
computed by Green's function method or other numeri- 
cal realization. The reason is that the numerical differ- 
ence, used for computing the current J = V x B from B, 
gives small but finite currents, whose directions are ran- 
domly from 0° to 180°, thus an average of the full volume 
should give a misalignment angle of ~ 90° . The noise in 
the observation data, mainly in the weak field regions, 
is also a major source for the random numerical errors. 
In these regions, the actual magnetic elements are proba- 
bly smaller than the observed or numerical pixel size, and 
the field directions generally exhibit a random pattern on 
the image. To reduce such errors in computing the met- 
ric, we can either put larger weight of current (e.g., use 
C^Wsin) or compute CWsin within the strong-current 
subregions only. As is expected, the misalignment an- 
gle decreases significantly by measuring in this way. For 
the full region, C^Wsin are only half of CWsin. For the 
subregions, CWsin are also only half or less for the full 
region, reaching the level of those from the benchmark 
tests (Jiang & Feng 2012a). In particular, the misalign- 
ment angle is only about 4° in subregion B, showing that 
the force-free assumption is modeled very well. 

Regarding the energy contents, it is interesting to note 
that the free energy of subregion A exceeds that of the 
full region, meaning that the free energy content in the 
full volume excluding subregion A is negative. This is, 



however, not surprising as we know that any sub- volume 
energy content of the non-potential field may be lower 
than the potential energy (e.g., Mackay et al. 2011; Jiang 
ct al. 2012a). Also the measurement error may result 
in this negative free energy since it is very small com- 
pared to the total free energy. No matter which case is 
true, it can be clearly seen that the spatial distribution of 
free energy is largely co-spatial with that of the current, 
since the subregion A contains most of the currents of 
the whole volume. This confirms that the free energy in 
the corona is actually stored by the current-carrying field 
(where non-potentiality is strong) , but not necessarily in 
the magnetic flux concentrations. 

Finally, we compare the results extrapolated from the 
raw and preprocessed magnetograms. Inspecting of the 
force- freeness and divergence- freeness metrics shows that 
the improvement by preprocessing is negligible. This 
is because the raw data already satisfies the boundary 
force-free conditions well. Due to the smoothing, the re- 
sult for the preprocessed data gives slightly lower energy 
contents than those for the raw data. 

4.2. AR 11283 

Figure 7 compares the AIA 171 A loops with the re- 
constructed field in the same way as Figure 5. For this 
AR, the NLFFF model appears to perform only slightly 
better than the potential model (there are some loops 
even worse produced by the NLFFF model than the po- 
tential model near the north-east boundary of the FoV). 
The clearest misalignment with the observation is the 
large closed loop pointed by the arrow in AIA image. 
This group of loops are failed to be recovered by both 
the potential and force-free models which give open field 
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Fig. 7. — Same as Figure 5 but for AR 11283. 
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Full 


0.32 
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5.12 


0.98 
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100% 
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0.13 


0.09 


1.76E-03 


2.05 


1.18 


0.87 
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TABLE 3 

Results of the metrics for AR 11283. Full region is extrapolation box of [0, 300](a;) x [0, 256](j/) x [0, 150](2). Region A is 

[193, 251] X [86, 140] x [0, 30]. The FoVs of the regions are shown by the box in Figure 8. The energy unit is 10'^^ erg. 
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Fig. 8. — Strong-current region and its magnetic structure for AR 11283. The upper-left panel is image of the vertical integral of current 
density, i.e., Jq = J \3\dz\ and upper-right panel shows the magnetic structure of subregion A, where a sigmoid, shown in the following 
panels, is observed clearly by SDO/AIA and Hinode/XRT. 
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lines instead. This, however, is not unexpected since the 
inputted magnetogram has flux unbalanced by —10%. 
So there must be field lines from the negative polarity 
opening in the FoV. The reason for this flux unbalance 
may be that the positive flux in the east is rather dis- 
persed (much more than the negative polarity), and thus 
properly be underestimated by the observation. 

Near the major polarity the structure of the loops is 
very complex and the extrapolated field shows highly 
sheared and twisted structures, indicating a significant 
non-potentiality there. Actually this was the site of the 
flare and filament eruption. The distribution of the ver- 
tically integrated current shows a strong concentration 
of current in this region, as denoted by A in Figure 8. In 
the same figure, we show the local field structure and the 
observations from different wavelengths of much higher 
temperature than 171 A. The magnetic field exhibits 
a multi-fiux rope configuration. The most remarkable 
structure is a sigmoid, i.e., the S-shaped loop in the AIA 
94 A and 335 A images. The sigmoid can be seen most 
clearly in the AIA 94 A wavelength (6.3 MK) with rather 
thin but enhanced shape, and is also well shaped in the 
soft X-ray image taken by Hinode/XKT. In the fourth 
panel of the figure the field lines are plotted overlying 
on the AIA 94 A image. It demonstrates that our ex- 
trapolation has recovered the sigmoid rather precisely, 
at least in morphology (see the precise alignment of the 
field lines with the shape of the sigmoid). The distribu- 
tion of the current also resembles roughly the shape of the 
sigmoid, suggesting that the enhancement of EUV and 
X-ray emission associated with the sigmoid is made by 
the strong field-aligned current via Joule heating of the 
plasma. This sigmoid locates between the major positive 
and negative polarities and the currents reside mostly in 
the north-east part, as shown by the current distribu- 
tion. The twist of the sigmoid field lines is not strong as 
modeled in other cases such as Rousscv et al. (2012) or 
Savcheva et al. (2012), and this sigmoid is composed of a 
single fiux rope, which is also different from their results 
with two flux ropes or double-J shaped current pattern. 
The observation and modeling suggest that there seems 
to be another flux rope overlying the sigmoid, and the 
flare and CME may be resulted by the eruptions of these 
flux ropes, which is left for future study. 

Similarly, we compute the metrics of force-freeness and 
divergence-freeness for both the full region and the sub- 
region and the results are given in Table 3. By comparing 
the results using the raw and preprocessed data, we flnd 
that evidently the preprocessed result is closer to force 
free, especially of the full region for which the raw data 
gives CWsin~ 0.4 (24°) while the preprocessed data gives 
CWsin~ 0.3 (17°). Thus for this AR the preprocessing 
indeed improves the extrapolation greatly. Also the di- 
vergence is reduced by the preprocessing. It is noticeable 
that the total energy content is doubled by the prepro- 
cessing, reaching 10'^^ erg. But even this improvement of 
the free energy is likely to underestimate the actual value, 
considering that a X-class flare and CME erupted imme- 
diately (Feng et al. 2013). Still the current is strongly 
localized and the free energy is concentrated within the 
strong-current region, i.e., subregion A, which occupies 
only less than one percent of the full volume, but contains 
most of the free energy. 



4.3. Convergence Study 

It is important to monitor the relaxation process to 
study whether the iteration converges, since there is no 
theory to guarantee this. Here we study the convergence 
process of the computations by temporal evolution of sev- 
eral monitors, including the residual of field between two 
successive iterations 



res"(B) = 



\ 



6—x,y,z 



E.iB-s)' 



(12) 



(where n denotes the iteration step), the metric CWsin 
and the total energy content. We record the residual by 
every ten steps and compute CWsin and total energy by 
every ten ta- Results for extrapolation of both ARs are 
plotted in Figure 9. As can be seen, the system con- 
verges smoothly and fast. During the first 10 ta, the 
residual keeps increasing because the transverse field is 
inputted at the bottom continuously, which drives the 
system away from the initial potential field. After this 
driving process, the residual drops immediately, indicat- 
ing a fast relaxation of the system. With about 40 ta 
(nearly 10000 iterations), the residual is already reduced 
to ~ 10^^, and all the metrics and energy almost stag- 
nate afterward. Thus the computations can actually be 
terminated once the residual is below 10~^, which is con- 
sistent with our previous studies for benchmark cases 
(Jiang & Feng 2012a; Jiang et al. 2012b). It is also note- 
worthy that the convergence process is rather smooth, 
without any obvious oscillation or abrupt variation of the 
residual or the metrics, so the iteration is "safe" . This is 
a good feature of our code over other iteration codes for 
extrapolation, e.g., the Valori et al. (2007) 's magnetofric- 
tional code or the Wheatland (2006) 's Grad-Rubin-like 
code, which usually show strong oscillatory in the itera- 
tion or even fail to converge occasionally (Schrijver et al. 
2008; DeRosaet al. 2009). 

5. CONCLUSIONS 

In this paper we have applied the CESE-MHD- 
NLFFF code to the SDO/BMI vector magnetograms. 
Two ARs are sampled for the test, AR 11158 and 
AR 11283, both of which produced X-class flares and 
were very non-potential. We compared the results with 
the SDO/AIA images, showing that the reconstructed 
field lines resemble well most of the plasma loops, which 
is a basic requirement for an applicable NLFFF model- 
ing code (DcRosa et al. 2009). Because the magnetic flux 
of the AR 11283 magnetogram is not well balanced, the 
extrapolation of the large scale fleld appears not as good 
as that for AR 11158. Observation shows that in the 
core regions of the ARs there were filament or sigmoid 
which are important precursors of eruptions like flares 
and CMEs. We also found in these places, there were 
highly-sheared and twisted field lines, i.e., flux ropes, 
which contain strong field-aligned currents and plenty 
of non-potential energy, and our extrapolations recov- 
ered indeed well those observed features, especially the 
sigmoid in AR 11283. By computing the metric CWsin 
which measures mean value of misalignment between the 
magnetic field and electric current, we found that, the 
force-free constraint is fulfilled very well in the strong- 
field regions (CWsin w 0.1, misalignment about 6°) but 



12 



Jiang et al. 



Iteration steps 
5000 10000 15000 



Iteration steps 
5000 10000 15000 



lO* 






lO* 






40 60 

t/T. 



t/T. 








1 1 1 1 


' ' 


1 


' 




\ 








- 


[ 
0.8 


B 

A 


Full 
A 








0.6 


:\ 








- 


0.4 


\^ 


3 B B 






- 


0.2 


: y 










- i 


^- i'-, ^ 


. 




i 






" 


" T 




, 


1 




1 


- 



20 40 

t/T. 






- 


A--=^=^ 


.--^ 


1 tu 1 


— UJ — 




T 


6 


i- ^-^ 






5.8 


- 


// 








- 


5.6 


- / 


7 








- 


5.4 


- / 


/ 










- 


















B— 


- Full 

- At4 




- 






bi^, 










1 


1 1 1 






1 


- 






20 


40 




60 







t/T, 



Fig. 9. — Convergence of computations: temporal evolution of the residual, the metric CWsin and the total energy content in the iteration 
process for AR 11158 (left column) and AR 11283 (right column). Note that for a better plot for the energy contents of different regions, 
their values are scaled properly as shown in the line legends. 
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apparently not that well in the weak-field regions (CWsin 
ft! 0.3, misalignment about 17°) because of the data noise 
and the numerical errors of the small currents. The en- 
ergy contents of our results are also consistent with the 
previous computations (with respect to the AR 11158, 
e.g., Wiegelmann et al. 2012; Sun ct al. 2012b). In sum- 
mary our extrapolation code can be used as a viable tool 
to study the 3D magnetic field in the corona. 

We developed the CESE-MHD-NLFFF code not only 
for field extrapolation, but also as a sub-program for 
the project of data-driven MHD modeling of the ARs, 
the eruptions and their dynamic evolutions in the global 
corona using continuously-observed data on the photo- 
sphere. At present numerical MHD investigations of the 
solar eruptions (Amari et al. 2003; MacNeice et al. 2004; 
Aulanicr ct al. 2009; Fan 2010; Torok ct al. 2011; Rousscv 
et al. 2012) are mostly based on idealized magnetic con- 
figurations without constrained by real observations. A 
step forward of understanding what really happens in the 
solar eruptions, certainly necessitates the observation- 
constrained numerical model. For example, considering 
that NLFFF extrapolation can recover highly-sheared 
magnetic arches and twisted flux ropes, which are basic 
building blocks of many eruption models (e.g., Torok & 



Kliem 2005; Aulanicr et al. 2009), utilizing the extrapo- 
lated field from real magnetograms can obviously provide 
much more realistic initial inputs than those idealized 
models like Titov & Dcmoulin (1999)'s flux rope model. 
Our future work is to input the extrapolation field as 
an initial condition into the data-driven full MHD model 
(Jiang et al. 2012a; Feng et al. 2012), along with the sur- 
face plasma flows derived from time-series of photosphere 
magnetograms (e.g., Liu et al. 2012) as bottom boundary 
condition to stress the model, with an objective to bet- 
ter simulate the initiation and evolution of solar explosive 
phenomena and their interplanetary evolution process. 
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